﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.EisPack
{
    /// <summary>
    /// Reduces a submatrix situated in rows and columnslow through high to upper hessenberg form 
    /// by stabilized elementary similarity transformations.
    /// </summary>
    [Serializable]
    public static class ElmhesClass
    {
        /// <summary>
        /// Reduces a submatrix situated in rows and columnslow through high to upper hessenberg form 
        /// by stabilized elementary similarity transformations.
        /// </summary>
        /// <param name="a">Contains the input matrix.</param>
        /// <param name="n">The order of the matrix.</param>
        /// <param name="lo">The lo parameter.</param>
        /// <param name="hi">The hi parameter.</param>
        /// <param name="intch">Contains information on the rows and columns interchanged in the reduction. 
        /// Only elements low through high are used.</param>
        public static void Elmhes(double[,] a, int n, int lo, int hi, int[] intch)
        {
            /* This routine is a translation of the Algol procedure from
             * Handbook for Automatic Computation, vol. II, Linear Algebra,
             * by Wilkinson and Reinsch, Springer-Verlag.
             */
            double x, y, dtmp;
            int i, j, kp1, la, m, mm1, mp1;

            la = hi - 1;
            kp1 = lo + 1;
            if (la < kp1)
            {
                goto _200;
            }
            for (m = kp1; m <= la; m++)
            {
                mm1 = m - 1;
                x = 0;
                i = m;
                for (j = m; j <= hi; j++)
                {
                    if (Math.Abs(a[j, mm1]) <= Math.Abs(x))
                    {
                        goto _100;
                    }
                    x = a[j, mm1];
                    i = j;
                    _100:
                    ;
                }
                intch[m] = i;
                if (i == m)
                {
                    goto _130;
                }

                /* Interchange rows and columns of array a. */
                for (j = mm1; j < n; j++)
                {
                    dtmp = a[i, j];
                    a[i, j] = a[m, j];
                    a[m, j] = dtmp;
                }
                for (j = 0; j <= hi; j++)
                {
                    dtmp = a[j, i];
                    a[j, i] = a[j, m];
                    a[j, m] = dtmp;
                }
                _130:
                if (x == 0.0)
                {
                    goto _180;
                }
                mp1 = m + 1;
                for (i = mp1; i <= hi; i++)
                {
                    y = a[i, mm1];
                    if (y == 0)
                    {
                        goto _160;
                    }
                    y /= x;
                    a[i, mm1] = y;
                    for (j = m; j < n; j++)
                    {
                        a[i, j] -= y * a[m, j];
                    }
                    for (j = 0; j <= hi; j++)
                    {
                        a[j, m] += y * a[j, i];
                    }
                    _160:
                    ;
                } /* endfor  i */
                _180:
                ;
            } /* endfor m */
            _200:
            ;
        }
    }
}